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ABSTRACT 

Simulation  of  orbital  phenomena  is  often  difficult  because  of  the  non-continuum  nature  of  the  flow, 
which  forces  the  use  of  particle  methods,  and  the  disparate  time  scales,  which  make  long  run  times 
necessary.  In  this  work,  the  computational  work  load  has  been  reduced  by  taking  advantage  of  the  low 
number  of  collisions  between  different  species.  This  allows  each  population  of  particles  to  be  brought 
into  convergence  separately  using  a  time  step  size  optimized  for  its  particular  motion.  The  converged 
populations  are  then  brought  together  to  simulate  low  probability  phenomena,  such  as  ionization  or 
excitation,  on  much  longer  time  scales.  The  result  of  this  technique  has  the  effect  of  reducing  run  times 

by  a  factor  of  1 0’  -  1 04.  The  technique  was  applied  to  the  simulation  of  a  low  earth  orbit  neutral  species 
sensor  with  an  ionizing  element. 


1.  Introduction 

The  measurement  of  neutral  species  levels  in  the  ionosphere  has  gained  considerable  importance 
due  the  effects  of  orbital  ionospheric  satellite  drag.  [1-4]  Variation  in  orbital  drag  can  change  satellite 
trajectories  and  result  in  difficulties  with  ground-to-orbit  communication,  In  addition,  satellite  motion 
and  thruster  firings  can  change  the  local  orbital  plasma  environment  leading  to  contamination  and 
alteration  of  fields  and  the  energies  of  particles  surrounding  the  spacecraft.[5-8]  Measurement  of  the 
neutral  component  of  the  ionosphere  is  usually  carried  out  with  some  type  of  mass  spectrometer.  These 
instruments  will  typically  first  ionize  the  incoming  flow  and  then  measure  the  current  to  characterize  the 
flux.  Once  the  particles  are  ionized,  electric  fields  can  then  be  used  to  extract  data  about  the  particle 
energy,  composition,  or  velocity  distribution  functioned]  In  an  effort  to  improve  the  design  of  an  ionizer 
element  for  a  neutral  particle  sensor,  numerical  simulation  can  be  used  to  optimize  the  performance  of  the 
instrument.  Simulation  ot  these  types  of  flow  fields  can  be  very  difficult  due  to  the  non-Maxwell  i  an 
nature  of  the  velocity  distribution  functions,  the  non-continuum  nature  of  the  fluid,  and  the  small 
ionization  fractions  involved.  All  of  these  problems  lead  to  large  computational  workloads  and  long 
computational  times.  The  present  work  details  the  strategies  and  techniques  utilized  in  the  acceleration 
ol  the  numerical  simulation  of  an  ionizing  element  used  in  an  orbital  neutral  particle  sensor.  The  purpose 
of  the  study  is  to  characterize  the  amount  of  current  which  can  he  expected  from  an  ionizing  element  used 
by  such  a  sensor. 

Computation  of  neutral  particle  flows  in  orbit  is  difficult  due  to  the  non-continuum  nature  of  the 
gas.  A  parameter  often  used  to  determine  the  expected  level  of  continuum  behavior  is  the  Knudsen 
number,  given  by, 


where  A  is  the  mean  free  path  and  L  is  the  length  scale  of  the  flow.  In  low  earth  orbit,  the  mean  free  path 
is  measured  in  thousands  of  meters. [10]  For  length  scales  comparable  to  a  satellite,  the  Knudsen  number, 
therefore,  is  greater  than  10'  and  the  flow  is  in  a  free-molecular  regime.  Therefore,  it  can  not  be  treated 
with  continuum  (fluid)  methods.  For  this  reason,  particle  methods  are  typically  employed  in  simulations 
ol  orbital  neutral  flows.  In  the  present  work,  a  Direct-Simulation-Monte-Cario[  1 1  ]  (DSMC)  technique  is 
utilized  to  model  particle  collisions  and  a  Particle-in-Ce!![12]  (PIC)  technique  is  used  to  model  the  effects 
of  electric  and  magnetic  fields. 

The  computations  will  attempt  to  simulate  the  ionization  process  in  an  orbital  neutral  particle 
sensor.  A  schematic  of  the  computational  geometry  is  shown  in  Fig.  1. 
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Fig.  1:  Computational  Geometry 
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The  neutrals  enter  the  computational  domain  on  the  left  side.  The  density  and  velocity  are  set  to  be 
representative  of  flux  of  particles  on  the  ram  side  of  a  satellite  which  is  in  an  orbit  with  a  600km  altitude. 
A  portion  of  the  neutrals  will  travel  through  a  gridded  grounding  screen.  For  the  present  simulations,  the 
screen  is  !  00%  transparent  and  serves  only  as  a  boundary  condition  for  the  electric  potential.  The  screen 
and  the  walls  of  the  ionization  chamber  are  grounded  to  the  spacecraft  chassis.  The  voltage  at  the 
detector  can  be  offset  from  ground  as  will  be  detailed  below.  Once  the  neutrals  enter  into  the  ionization 
chamber,  they  are  bombarded  by  a  beam  of  electrons  from  a  source  located  on  the  top  of  the  ionization 
chamber.  The  energy  of  the  electrons  is  set  to  lOOeV  to  match  the  peak  of  the  ionization  cross  section  for 
atomic  oxygen,  the  dominant  neutral  in  LEO.  The  resulting  ions  are  then  collected  in  the  detector  located 
on  the  right  side  of  the  ionization  chamber.  The  computations  treat  the  geometry  in  a  2-dimensional 
manner,  although,  the  particles  have  3-dimensional  velocity  distributions.  The  geometry  is  assumed  to  be 


!cm  in  the  direction  normal  to  the  page.  The  purpose  of  the  study  is  to  determine  the  amount  of  ions 
which  can  be  expected  to  reach  the  sensor  after  being  ionized  by  the  electron  beam. 


2.  Method 

The  simulation  of  the  flow  field  detailed  above  is  straight  forward  but  is  computationally  very 
intensive.  It  requires  modeling  of  three  major  physical  phenomena:  I)  the  neutral  flow  and  associated 
collisions,  2)  the  electron  flow,  3)  the  ionizing  collisions,  and  4)  the  motion  of  the  created  ions.  By 
inspecting  each  of  these  phenomena  and  dispensing  with  the  components  which  do  not  affect  the  results 
in  a  statistically  significant  manner,  the  computational  workload  (run  time)  can  be  reduced  bv  orders  of 
magnitude. 

We  will  first  examine  the  neutral  flow  and  associated  collisions.  The  neutral  flow  will  enter  with 
a  non- Maxwellian  velocity  distribution.  For  the  purposes  of  this  study,  the  neutral  velocity  distribution 
m  the  X  (ram)  and  the  Y  and  Z  (perpendicular)  directions  wili  be  as  shown  in  Fig.  2.  The  X  velocity 
distribution  represents  an  orbital  velocity  flow  of  7600m/s  with  a  temperature  of  2000K  while  the  Y  and 
Z  velocity  distributions  represent  a  neutral  population  with  a  temperature  of  1 000K.  Each  distribution  is 
essentially  a  ID  drifted  Maxwellian,  although  the  Y  and  Z  components  have  no  directed  velocity.  The 
different  temperatures  are  chosen  to  be  representative  of  a  flow  which  may  have  been  accelerated  in  an 
external  electric  field.  In  the  computations,  these  velocity  distribution  functions  (VDF)  are  sampled  in 
order  to  generate  the  population  of  neutrals  entering  the  computational  domain. 


After  entering  the  computational  domain,  the  neutral  particles  can  collide  with  other  neutral  particles. 

I  he  average  neutral  collision  frequency  can  be  estimated  with, 

Ko,=nn(rVlh  (2) 

Where  n„  is  the  neutral  density,  a  is  the  collision  cross  section,  and  V„,  is  the  thermal  velocity.  For  the 
current  conditions,  the  neutral-neutral  collision  frequency,  vL.oh  is  0.015  s'1.  The  transit  time  for  the 


neutrals  through  the  computational  domain  (moving  at  orbital  velocity  of  7.6km/s)  is  only  about  IfT5  s. 
Therefore,  the  probability  that  a  neutral  collides  with  another  neutral  is  of  the  order  of  10"7  or  only  one  in 
10  million  neutrals  will  suffer  a  collision.  If  the  neutral-neutral  collisions  were  to  be  ignored,  it  would 
only  incur  an  error  in  the  velocity  distribution  function  of  one  part  in  1 07.  For  this  reason,  neutral-neutral 
collisions  will  not  be  included  in  the  computations.  This  will  save  the  computational  cost  of  calculating 
probabilities  and  pairing  particles  which  must  still  be  carried  out  even  though  no  significant  number  of 
neutral-neutral  collisions  will  take  place.  Neutral  particles  will  still  be  allowed  to  collide  with  walls  and 
suffer  ionizing  collisions  with  electrons  but  neutral-neutral  collisions  will  not  be  calculated.  Effectively, 
their  velocity  distribution  functions  will  only  be  changed  by  wall  collisions. 

I  he  second  phenomena  to  be  considered  is  the  electron  flow.  This  includes  electron  transport, 
electron-electron  collisions,  momentum  transfer  electron-neutral  collisions,  and  electron-neutral  ionizing 
collisions.  Electron  transport,  the  acceleration/deceleration  and  movement  of  electrons  in  ambient  electric 
fields,  is  an  important  first  order  effect  and  will  be  retained  in  the  computations.  Electron-electron 
collisions,  however,  do  not  affect  the  electron  velocity  distribution  function.  This  is  because  two 
colliding  electrons  have  identical  mass  and  the  result  of  conserving  momentum  and  energy  is  that  the 
electrons  simply  exchange  their  velocity  vectors.  Therefore,  there  is  no  net  effect  on  the  electron 
population  velocity  distribution.  For  this  reason,  electron-electron  collisions  also  will  not  be  included  in 
the  computations. 

Next,  electron-neutral  momentum  transfer  collisions  will  be  considered.  The  electrons  will  enter 
the  simulation  at  a  rate  of  0. 1  mA  with  an  energy  of  lOOeV  and  a  temperature  of  1 000K.  The  average 
energy  is  chosen  to  match  the  peak  of  the  atomic  Oxygen  electron  impact  ionization  cross-section.  At  this 
energy,  the  average  electron  will  be  travelling  with  a  velocity  of  5.9x1 06  m/s  and  will  cross  the  ionization 
chamber  (2.5cm  tall)  in  4x  1 0  s.  Conservatively  estimating  the  average  elastic  electron-neutral  cross 
section  to  be  of  the  order  of  1 0  l9m2  and  using  eqn.  2  gives  a  collision  frequency  of  about  6  per  second. 
Since  the  electrons  will  only  be  in  flight  for  1G"9  s,  the  probability  of  an  electron  suffering  an  elastic 
collision  with  a  neutral  particle  is  10  s,  or  one  collision  for  every  1 00  million  electrons.  For  this  reason, 
the  computations  will  also  ignore  electron-neutral  momentum  transfer  collisions. 

Next,  electron  impact  ionization  collisions  will  be  considered,  The  atomic  Oxygen  electron 
impact  ionization  cross-section  at  lOOeV,  the  average  electron  energy,  is  about  1 .4x1  O'20  m2.  Again,  using 
eqn.  2,  this  results  in  an  average  ionization  collision  frequency  of  0.83s'1.  When  coupled  with  the 
electron  transit  time,  the  average  probability  of  ionization  is  1 0'9.  This  presents  a  tremendous 
computational  problem:  on  average,  109  electrons  must  be  simulated  for  every  ionization  event  that  takes 
place.  This  is  exceedingly  small.  However,  this  phenomena  can  not  be  ignored  since  determining  the 
amount  of  ionization  is  the  principal  goal  of  the  study. 

The  present  ionizer  will  use  an  electron  current  of  0. 1  mA.  This  represents  1014  electrons  per 
second  which  would  produce  1 06  ions/sec  using  the  above  stated  probability.  In  order  to  take  statistics 
on  I04  ions,  a  total  simulation  time  of  1  O'2  seconds  is  needed.  Since  the  electrons  are  travelling  at  lOWs 
and  flic  ceU  size  is  of  the  order  of  I0'3m,  the  electrons  will  be  limited  to  time  step  of  the  order  of  I0'9- 
10  s.  Therefore,  in  order  to  simulate  a  total  of  10  "s,  a  total  of  1 0s  time  steps  are  needed.  This  represents 
an  unacceptable  computational  work  load.  Parallelization  of  the  computer  code  will  not  alleviate  this 


problem  since  the  flow  will  be  evolving  in  time  and  must  be  computed  in  chronological  order.  Clearly,  a 
different  computational  strategy  is  needed. 


We  can  take  advantage  of  the  physics  of  the  current  problem  in  order  to  decrease  the  work  load 
by  noting  the  influence  that  each  phenomena  has  on  the  others.  As  explained  above,  because  electron- 
neutral  momentum  transfer  collisions  are  so  infrequent,  the  neutral  particle  velocity  distribution  function 
is  largely  independent  of  the  electron  population  by  a  factor  of  better  than  one  part  in  10s.  This  suggests 
that  the  neutral  population  in  the  computational  domain  can  be  simulated  independent  of  the  electrons. 
This  would  allow  the  neutrals  to  be  brought  to  equilibrium  using  a  time  step  that  is  consistent  with  the 
neutrals  and  not  restricted  by  the  high  velocity  of  the  electrons.  For  the  current  conditions,  the  neutral 
time  step  is  determined  by  the  orbital  velocity  (7600m/s)  and  the  ceil  size  (KTm)  and  is  of  the  order  of 
1 0'7s,  a  1 000  fold  speed  up  over  the  simulation  which  includes  electrons.  Similarly,  since  the  electron 
distribution  is  independent  of  the  neutrals,  their  equilibrium  condition  can  be  computed  without  having  to 
compute  neutral  motion.  After  independent  convergence  to  equilibrium,  the  neutral  and  electron 
populations,  with  the  imbedded  velocity  distribution  information,  can  be  saved  and  used  as  input  to  a 
simulation  involving  only  the  ionization  and  transport  of  the  ions.  If  the  ion  population  remains  small  in 
comparison  to  the  electron  population,  the  interaction  between  the  two  need  not  be  computed.  This  final 
simulation  would,  therefore,  not  need  to  compute  any  changes  to  the  neutral  or  electron  populations.  This 
can  not  be  determined  a  priori,  but  will  turn  out  to  be  the  case  as  will  be  shown  below. 


In  order  to  improve  statistics  and  further  reduce  run  times,  one  further  modification  is  proposed. 
Since  the  time  step  will  need  to  be  as  small  as  10'7s  in  order  to  follow  the  ion  motion,  an  ionization  rate  of 
1 0  ’  per  second  will  not  produce  many  ions  during  each  time  step.  Since  statistical  significance  increases 
with  the  number  of  computational  particles  in  the  simulation,  more  ion  particle  creation  during  each  time 
step  is  desired.  This  can  be  achieved  by  allowing  the  creation  of  fractional  ions  or  computational  particles 
that  represent  only  a  fraction  of  an  ion.  This  will  be  accomplished  in  the  following  manner:  The 
precomputed  neutral  and  electron  populations  will  be  used  to  determine  the  total  ionization  rate  in  each 
cell.  During  each  time  step,  a  single  particle  representing  a  fraction  of  an  ion  will  be  generated  in  each 
cell.  The  velocity  of  the  ion  will  be  determined  by  carrying  out  a  collision  between  a  randomly  chosen 
neutral  and  electron  particle.  The  weighting  of  the  particle,  however,  will  be  set  to  the  number  (or 
fraction)  ot  ions  which  statistically  should  be  generated  in  each  cell  during  each  time  step.  Therefore, 
instead  of  waiting  an  average  of  10-  1 0,000  time  steps  for  each  new  ion,  ions  will  be  generated  during 
each  and  every  time  step  and  their  population  and  behavior  will  be  determined  in  a  fractionally  weighted 
manner.  Fractional  particles  will  not  change  the  physics  since  the  density  will  be  low  enough  to  preclude 
significant  numbers  of  collisions.  In  effect,  the  scheme  remains  unchanged  from  where  a  single 
computational  particle  represented  multiple  real  particles.  This  strategy  will  generate  more  computational 
ion  particles  but  the  same  average  number  of  real  physical  particles.  Therefore  the  statistics  will  be 
improved  and  the  required  run  times  will  be  reduced. 


The  effect  on  the  electron  population  still  needs  to  be  determined.  If  the  ion  density  magnitude 
approaches  the  electron  density  magnitude,  the  electric  fields  will  begin  to  be  cancelled  and  the  electron 
flow  will  need  to  be  recomputed.  The  ion  density  will  be  increased  by  the  ionization  events  and 
decreased  by  the  neutralization  at  the  walls.  The  ion  flow  velocity  will  start  near  the  neutral  orbital 
velocity  (because  the  ions  are  born  with  the  velocity  of  the  parent  neutral)  and  then  may  change  as  the 


ions  are  accelerated  by  local  fields.  Therefore,  there  is  no  way  to  determine  the  final  ion  density  a  priori. 
Computations,  which  will  be  shown  below,  indicate  that  the  ion  density  never  rises  above  107  -  I0®m'3. 
Comparing  this  to  the  average  electron  density  (10,2m'3)  shows  that  there  will  exist  104  -  10s  electrons  for 
every  ion.  1  herefore,  they  will  not  significantly  change  the  potential  distribution  inside  the  ionization 
chamber  and  the  effect  on  the  electrons  will  not  be  significant.  In  addition,  the  low  density  of  the  ions 
results  in  a  very  small  collision  frequency  between  them  and  the  neutrals  (  vC0|  =  0.006  s'').  During  the 
transport  to  the  walls  (about  1 0'6s  at  near  orbital  velocity),  the  probability  of  an  ion  suffering  a  collision  is 
of  the  order  of  10'9.  Therefore,  since  only  one  particle  in  \09  will  sufFer  a  collision,  the  ion  flow  will  not 
impact  the  neutral  distribution  either  and  can  be  computed  independently  of  both  the  electron  and  neutral 
particle  distributions. 

In  summary,  because  of  the  low  collisionality  between  neutrals  and  electrons,  each  will  be 
computed  separately  and  using  a  time  step  optimized  to  its  motion,  The  converged  populations  will  then 
be  frozen  and  used  in  the  ionization  source  terms.  This  will  allow  computation  of  the  low-probability 
ionization  events  over  long  time  periods  without  the  necessity  of  computing  the  neutral  and  electron 
motion.  Because  of  the  low  collisionality  between  the  ions  and  the  neutrals,  the  resulting  ion  flow  will 
also  be  computed  separately.  In  order  to  improve  statistics  and  further  reduce  run  times,  fractional  ion 
generation  is  employed  which  reproduces  the  statistical  ionization  rate.  The  result  of  these  techniques  is  a 
reduction  in  the  number  of  time  steps  necessary  to  reach  convergence  from  1 0s  to  a  number  less  than  10s. 
This  has  the  effect  of  reducing  the  run  times  by  a  factor  of  1 03  -  104. 


3.  Results 


The  computations  are  initiated  using  the  geometry  detailed  in  Fig.  1 .  The  computational  grid  is 
formed  by  cells  which  are  1mm  long  (X)  by  I  mm  tall  (Y)  by  10mm  into  the  page  (Z).  This  represents  a 
cell  size  many  orders  of  magnitude  smaller  than  the  mean  free  path  (>]  03  m)  and  smaller  than  the  Debye 
length  (—  0.004  m).  Note  that  the  Debye  length  is  strictly  not  applicable  in  the  current  simulations 
because  the  flow  is  not  a  true  plasma  since  it  is  never  quasi-neutral.  Therefore,  it  never  exhibits  the 
shielding  behavior  ol  an  actual  plasma.  Never-the-less,  the  cel!  size  is  kept  small  in  order  to  maintain 
sufficient  accuracy  in  modeling  of  the  geometry.  The  computations  arc  treated  as  two  dimensional 
although  the  particles  have  velocities  in  the  three  dimensions.  Particles  leaving  the  domain  in  the  Z 
direction  (into  the  page)  simply  reenter  the  domain  from  the  opposite  side  maintaining  the  same  Z 
velocity.  Particles  striking  walls  are  either  reflected  in  a  diffuse  or  specular  manner.  With  diffuse 
reflection,  the  particles  striking  the  wall  leave  with  direction  to  the  wall  that  has  a  cosine  probability- 
distribution.  The  average  velocity  is  determined  by  the  wall  temperature  (full  accommodation)  which  is 
set  to  300  K.  With,  specular  reflection,  the  particles  leave  the  wail  with  the  same  angle  to  the  normal  as 
they  had  on  impact.  Ail  neutrals  enter  the  simulation  from  the  left  side  of  the  domain  with  a  velocity 
distribution  given  by  Fig.  2.  The  velocity  distribution  is  further  modified  such  that  the  x  and  y  velocities 
are  rotated  1  Odeg.  This  simulates  a  sensor  that  is  pointed  in  a  direction  that  is  not  exactly  into  the  RAM. 
Computations  are  continued  for  several  hundred  neutral  flight  times  which  allows  the  number  of  particles 


in  the  domain  to  stabilize.  The  resulting  neutral  density  contours  for  diffuse  and  specular  wall  reflection 
are  displayed  in  Fig.  3. 


a)  diffuse  wall  reflection  b)  specular  wall  reflection 


Fig.  j.  Neutral  number  density  contours  (m J)  tor  diffuse  and  specular  wall  reflection  cases.  The  ambient 
entry  density  is  lxlO'mv’  with  the  flow  entering  at  an  angle lOdeg.  upward. 

Computations  show  that  the  result  of  the  neutrals  entering  with  a  velocity  angled  upward  by  10 
deg.  is  a  neutral  pile  up  near  the  top  of  the  ionization  chamber,  especially  for  the  specular  reflection  case. 
Hie  computations  also  show  that  diffuse  reflection  results  in  higher  density  levels  in  the  ionization 
chamber  of  the  sensor.  The  average  density  for  the  diffuse  reflection  case  is  40  times  the  ambient  density, 
whereas  for  specular  reflection,  it  is  only  a  factor  of  2  times  the  ambient  density.  This  occurs  because, 
with  diffuse  reflection  and  full  accommodation,  the  particles  leave  the  wall  with  a  much  smaller  velocity 
(wall  thermal  vs  orbital  or  7Q0m/s  vs  7600m/s)  since  their  energy  is  set  to  equal  the  wall  temperature. 

Ref.  1 3  suggests  that  at  LEO,  diffuse  reflection  may  be  a  more  appropriate  model,  while  at  higher 
altitudes  such  as  GEO,  specular  reflection  may  dominate.  In  order  to  make  a  conservative  estimate  of  the 
lowest  number  of  ions  impacting  the  detector,  this  study  will  use  the  specular  reflection  model  and  the 
resulting  lower  levels  of  neutral  density. 


The  electron  flow  was  simulated  independently  of  the  neutral  flow  because  of  the  low  number  of 
collisions  taking  place  between  the  two,  as  detailed  in  the  previous  section.  The  electric  potential  was 
solved  at  each  time  step  and  the  electric  field  forces  were  applied  to  the  electron  particles  as  they  moved 
through  the  domain.  Ail  walls  were  held  at  a  0  voltage  boundary  condition.  Electrons  impacting  the 
walls  were  removed  from  the  computation  because  the  walls  were  modeled  as  being  made  of  a  conducting 
material.  No  secondary  electron  emission  is  currently  modeled.  Three  electron  current  levels  were 
simulated:  0.01,  0. 1 ,  and  1  mA.  The  resulting  electron  number  density  contours  are  show  in  Fig.  4. 


a)  Ic  =  G.OlrnA 


b)lc“  0.1mA  c)Ie=  1.0mA 


1"  ig.  4,  Electron  number  density  contours  (m  '  )  for  different  electron  current  levels. 


As  the  electron  current  increases  from  0.01  -  1 .0  mA,  the  average  electron  number  density 
increases  linearly  from  I0n  -  10"  m"J.  The  electron  beam  width  also  increases  with  increasing  current 
levels  and  in  the  case  where  the  current  reaches  1 .0  mA,  the  electron  beam  clearly  is  impinging  on  the 
detector.  This  anomalous  electron  flux  may  alter  the  measured  current  at  the  detector.  Current 
contamination  could  be  avoided  by  simply  moving  the  electron  emitter  to  a  position  farther  from  the 
detector  or  biasing  the  detector  face.  The  reason  for  the  beam  broadening  can  be  found  in  the  electric 
potential  contours  shown  in  Fig.  5. 


Fig.  5  shows  that  the  presence  of  the  electron  beam  creates  a  negative  potential  region  in  the 
ionization  chamber.  The  level  of  negative  potential  increases  as  the  electron  current  increases,  from  about 
-0.05  V  for  a  current  of  0.0 1  mA  to  -5  V  at  an  electron  current  level  of  1 ,0mA.  The  negative  potential  well 
creates  electric  fields  that  spread  the  electron  beam.  This  negative  potential  well  also  creates  another 
problem:  ions  born  inside  it  must  have  sufficient  energy  to  climb  out  and  reach  the  detector.  For  the 


present  case,  the  atomic  oxygen  neutrals  enter  with  an  energy  of  4.6eV  due  to  the  7600m/s  entry  speed. 
Since  the  bulk  of  the  energy  of  ionization  will  be  carried  away  by  the  liberated  electron,  ions  will  be  born 
with  an  energy  close  to  that  of  the  neutrals.  An  ion  with  less  than  5eV  will  not  be  able  to  reach  the 
detector  it  it  is  bom  at  the  center  of  the  potential  well  created  by  the  presence  of  the  electron  beam. 

I  he  neutral  and  electron  particle  populations  computed  above  were  then  used  to  calculate 
ionization  rates  in  a  simulation  of  the  ions.  Neutral  and  electron  particles  were  chosen  at  random  in  each 
cell  in  order  to  preserve  the  collisional  dynamics  of  the  electron  impact  collisions  and  the  resulting  ions 
were  allowed  to  move  under  the  influence  of  the  local  electric  fields.  The  resulting  ion  number  density 
contours  are  shown  in  Fig.  6. 


Fig.  6:  Ion  number  density  contours  (m  ’)  for  different  electron  current  levels. 


I  he  results  show  that  the  average  ion  density  increases,  as  expected,  with  increasing  electron 
current.  In  each  case,  ions  are  noted  to  exit  the  ionizing  chamber.  This  flux  ofions  is  due  to  ions  created 
from  neutrals  that  had  already  rebounded  from  the  detector  face  and  had  a  velocity  pointed  outward  from 
the  detector.  For  the  two  lowest  electron  current  levels,  the  ion  density  increases  linearly  with  electron 
current.  For  the  last  case  (Ie  =  1.0mA),  however,  the  ion  density  increases  more  than  can  be  attributed  to 
the  increased  electron  current  alone.  In  this  case,  the  negative  potential  well  (Fig.  5c)  is  of  the  same  order 
as  the  neutral  energy.  This  results  in  a  dynamic  situation  where  a  large  portion  of  the  ions  inside  the 
ionization  chamber  are  trapped  due  to  their  inability  to  climb  out  of  the  potential  well.  These  ions  would 
presumably  be  released  and  be  free  to  impact  the  detector  or  leave  the  ionizing  chamber  once  the  electron 
current  was  turned  off  and  the  resultant  potential  well  no  longer  existed.  Therefore,  operating  the  sensor 
at  high  electron  current  levels  can  create  situations  where  a  significant  portion  of  the  created  ions  can  not 
reach  the  detector. 


The  problem  of  ion  trapping  can  be  mitigated  by  biasing  the  detector  face.  For  example,  if  the 
detector  in  the  1 ,0mA  case  were  to  be  biased  to  a  level  higher  than  the  value  of  the  potential  well,  a  path 
would  be  created  where  the  ions  can  escape  and  reach  the  detector.  This  was  attempted  by  biasing  the 
detector  face  to  -1 5  V  and  repeating  the  computations.  The  resulting  electric  fields  did  not  affect  the 


neutral  flow,  therefore,  it  was  unnecessary  to  recalculate  them.  The  bias  was  applied  and  the  electrons 
and  ions  were  recalculated.  The  results  are  shown  in  Fig.  7. 


Fig.  7:  Effects  on  electric  potential,  electron  density,  and  ion  density  of  biasing  the  detector  face. 


Figure  7a  shows  the  electric  potential  inside  the  ionization  chamber  when  the  detector  face  is 
biased  to  -15V  and  a  current  of  ImA  is  used.  There  is  now  a  zero  voltage  potential  difference  between 
the  center  of  the  ionization  region  and  the  detector.  This  allows  all  of  the  created  ions  to  reach  the 
detector.  1  he  electron  density  distribution,  shown  in  Fig.  7b,  shows  more  spreading  than  exhibited  in  Fig. 
4c.  This  is  due  to  the  electric  fields  (potential  gradient)  created  by  the  biased  detector.  The  average  ion 
density,  shown  in  fig.  7c,  now  reaches  only  10'm  ’  which  continues  the  linear  variation  with  electron 
current  exhibited  by  the  0.01  and  0.1mA  cases.  This  indicates  that  ions  are  no  longer  being  trapped.  In 
fact,  the  lack  of  ions  outside  of  the  ionization  chamber  in  Fig.  7c,  suggest  that  even  the  ions  bom  with 
velocities  which  would  normally  take  them  out  of  the  chamber  are  now  largely  being  attracted  to  the 
detector.  1  he  ion  impacts  to  the  detector  are  shown  in  Fig.  8.  At  low  electron  currents,  0.01  and  0. 1  mA. 
the  potential  well  created  is  smaller  than  the  ion  energy  and  biasing  the  detector  does  not  affect  the  ion 
current  to  the  detector.  At  an  electron  current  of  ImA,  however,  biasing  the  detector  face  is  necessary  to 
collect  all  of  the  created  ions  and  avoid  ion  trapping. 


detector 


ettctron  current  (mA) 

Fig*  8:  Effects  of  detector  bias  on  the  number  of  ion  impacts  at  the  detector. 


4*  Conclusions 


Simulation  of  orbital  phenomena  is  often  difficult  because  of  the  non-continuum  and  non-Maxwellian 
nature  of  the  flow,  which  forces  the  use  of  particle  methods,  and  the  disparate  time  scales  from  low 
probability  phenomena,  which  make  long  run  times  necessary.  In  this  work,  the  computational  work  load 
has  been  reduced  by  taking  advantage  of  the  low  number  of  collisions  between  different  species.  This 
allow  s  each  population  of  particles  to  be  brought  into  convergence  separately  using  a  time  step  size 
optimized  for  its  particular  motion.  The  converged  populations  can  then  be  brought  together  to  simulate 
low  probability  phenomena  on  much  longer  time  scales.  The  result  of  this  technique  has  the  effect  of 
reducing  run  times  by  a  factor  of  I03  -  104.  The  technique  was  applied  to  the  simulation  of  a  low  earth 
orbit  neutral  species  sensor  with  an  ionizing  element.  The  simulations  demonstrate  the  viability  of  such  a 
sensor  and  reveal  the  need  for  detector  face  bias  when  large  ionizing  currents  are  used. 


References 

[1]  M,  F.  Storz,  B.  R.  Bowman,  J.  i.  Branson,  S.  J.  Casali,  and  W.  K.  Tobiska,  High  Accuracy  Satellite 
Drag  Model  (HASDM),  Advances  in  Space  Research,  36(12)  (2005)  2497-2505. 

[2]  J.  F.  Liu,  R.  G.  France,  and  H.  B.  Wackemagel,  An  Analysis  of  the  Use  of  Empirical  Atmospheric 
Density  Models  in  Orbital  Mechanics,  AAS/AIAA  Astrodynamics  Specialist  Conference,  (1983). 

[3]  F.  A.  Marcos,  Accuracy  of  Atmosphere  Drag  Models  at  Low  Satellite  Altitudes,  Adv.  Space  Research. 
10(3)  (1990)  4!  7-422. 

[4]  Nicholas,  I .  Finne,  I.  Galysh,  M.  Davis,  and  L.  Healy,  The  Atmospheric  Neutral  Density  Experiment 
(ANDE),  2009  NRL  Review  (2009)  138-142. 


[5]  E.  Wulf  and  U.  vonZahn  .  The  shuttle  environment:  Efreets  of  thruster  firings  on  gas  density  and 
composition  in  the  payload  bay,  J.  Geophys,  Res.,  91{A3),  (1986)  3270-3278, 
doi:  1 0. 1 029/J  AQ9 1  iA03n03270. 

[6j  J.  S.  Machuzak,  W.  J.  Burke,  J.  M.  Retterer,  D.  E.  Hunton,  J.  R,  Jasperse,  and  M.  Smiddy,  Effects  of 
thruster  firings  on  the  shuttle’s  plasma  and  electric  field  environment,  J.  Geophys.  Res.,  98(A2)  (1993) 

1 5 1 3-1 530,  doi:  1 0. 1 029/92 JA0204 1 . 

[7]  U.  Samir,  N.  H.  Stone,  K.  H.  Wright  Jr.,  On  plasma  disturbances  caused  by  the  motion  of  the  space 
shuttle  and  small  satellites:  A  comparison  of  in  situ  observations,  J.  Geophys.  Res.,  91  (1986)  277-285. 

[8]  U.  Samir,  and  G.  L,  Wrenn,  Expermental  Evidence  of  an  Electron  Temperature  Enhancement  in  the 
Wake  of  an  ionospheric  Satellite,  Planet.  Space  Sci.,  20  (1972)  899-904. 

[9]  C.  Enloe,  K.  Habush,  R.  Haaland,  T.  Patterson,  C.  Richardson,  C.  Lazidis,  and  R.  Whiting, 
Miniaturized  electrostatic  analyzer  manufactured  using  photolithographic  etching,  Rev.  Sci.  Instrum., 
74(3) (2003) 1192-1195. 

1 1 0 j  D.  Hastings,  and  H.  Garrett,  Spacecraft-Environment  Interactions,  Cambridge  University  Press, 
Cambridge,  1996, 

[HI  G.  Bird,  Molecular  Gas  Dynamics  and  the  Direct  Simulation  of  Gas  Flows,  Oxford  University  Press, 
Oxford,  1994. 

112]  C.  Birdsall,  and  A.  Langdon  ,  Plasma  Physics  via  Computer  Simulation,  Tavlor  and  Francis,  New 
York,  2004. 

[13]  K.  Moe,  and  M.  M.  Moe.  Gas-surface  interactions  and  satellite  drag  coefficients.  Planet.  Space  Sci.. 
53  (2005)  793-801.  DOI:  10. 1016/j.pss.2005.03.005 


